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Abstract 

Soft-core attractive potentials can give rise to a phase diagram with three fluid phases at 
different densities (gas, low-density liquid and high-density liquid), separated by first order phase 
transition lines ending in critical points. Experiments show a phase diagram with these features 
for phosphorous and triphenyl-phosphite. Liquid-liquid phase transition could be relevant for 
water, silica, liquid metals, colloids and protein solutions, among others. Here we compare two 
potentials with short-range soft-core repulsion and narrow attraction. One of them is a squared 
potential that is known to have liquid-liquid phase transition, ending in a critical point, and no 
anomaly in density. The normal, monotonic, behavior of density for isobaric cooling is surprising 
if compared with molecular liquids, such as water, where a hypothetical critical point is proposed 
as rationale for the anomalous behavior of density. The second potential is a continuous version 
of the first. We show that the phase diagram associated to this new potential has, not only the 
liquid-liquid phase transition, but also the density anomaly. Our result, therefore, shows that the 
behavior in density is strongly dependent on the derivative of the potential. 
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I. INTRODUCTION 



Increasing experimental, theoretical and numerical evidences are showing that single- 
component systems can have more than the two commonly known fluid phases: gas and 
liquid. There is no thermodynamic inconsistency, indeed, in hypnotizing the possible ex- 
istence of a phase transition between two (or more) liquids at different densities. Direct 
experimental evidences of this phenomenon have been observed in 
and triphenyl phosphite 3], while consistent data exists for water [4 



aluminate liquids 



111 ] , selenium [12| , and cobalt |l3|, am ong others 



a liquid-liquid critical point for supercooled water 



221 ] , supercooled silica |9j, l23l. 



16L 



literature is contradictory 



26. 



27 



and hydrogen 



15 



16 



18 



iquid phosphorous m \2 



6j, iTJ, |8| , silica 



14J . Simulations predict 



20 



21], phosphorous 



251 ] . For other substances, such as carbon, 



28]. 



In the last years, several models have been proposed to understand the origin of the 
liquid-liquid phase transition. Within this context a large interest has been focused on 
isotropic core-softened models [291]. The reason is twofold. On the one hand, within accept- 
able limitations, they are models for a variety of systems including liquid metals, metallic 



mixtures, electrolytes, colloids and 



and silica 



30 



31 



33 



34 



35 



protein so 



36 



37 



38 



ution, as well as anomalous liquids, like water 



39 



40 



4 II ]. On the other hand, their sim- 



ple definition encourages to investigate the possibility of an intriguing relation between the 
liquid-liquid phase transition and the anomalies in specific observables, such as the density 
or the diffusivity, in network-forming liquids, like water and silica. 

Different kinds of isotropic potentials with soft-core have been proposed. In many of them, 
with the soft-core given by a Gaussian core 4^| or a ramp, the presence of a liquid-liquid 
phase transition is observed together with anomalies in the density, in the diffusion and in the 



structure 



39 



40j, |43j, M, 



45 



46 



47J, |48j. This is consistent with what is proposed 



on the basis of simulations and theoretical calculations for specific models 



m 



18 



or water 
21, 291 



In other attractive isotropic potentials, with the soft-core given by a square s 



the liquid-liquid phase transition occurs in absence of density anomaly 49 



50 



roulder, 



51 



52]. 



Nevertheless, an asymptotic tendency to the density anomaly is observed in three dimensions 



3 



and thermodynamic and dynamic anoma 



dimensions on lattices 



53j or off- lattice |45|,|54|. 



ies has been found in related models in two 



Here we focus on the study of a potential that is the smooth version of the one with 



2 



square shoulder and attractive well which we studied in Ref. 49j. We find that the new 
continuous potential has a liquid-liquid phase transition, as well as the previous squared 
potential, consistent with the minimal difference between the two models. Nevertheless, the 
difference is enough to cause the appearance of the density anomaly. Hence, the anomalous 
behavior is strongly dependent on the fine details of the potential, including its derivative. 



II. THE CONTINUOUS POTENTIAL 

We consider the following isotropic pairwise potential 



where Ur and Ua are the energy of the repulsive shoulder and of the attractive well, re- 
spectively; a, Rr and Ra are the hard-core distance, the repulsive average radius and the 
distance of the attractive minimum, respectively; A is a parameter related to the slope of 
the potential at R R ; 5 A is the variance of the Gaussian centered at Ra- 

For the set of parameters Ur/Ua = 2, Rr/cl = 1.6 Ra/o, = 2, A = 15, = 0.1, 

at distance r/a = 3 the value of the potential is Uq/Ua — —6.7 x 1CT 3 , with Uq — > for 
r/a > 3. Hence, to reduce the computational effort, we set a potential cutoff at r/a = 3 and 
add a constant and a linear term to the potential in order to have both the potential and 
its r-derivative equal to zero at the cutoff. 

The resulting potential is 

U(r) = U (r) + CU A + B- (2) 

(X 

where, for our choice of parameters and cutoff distance r c /a = 3, the dimensionless constants 
are set to C — 0.208876 and B = —0.0673794. The approximate attractive width, calculated 
as the width at half height of the Gaussian function defining the attractive well, is wa = 
25 AV /21n(2) ~ 0.75 a. 

The continuous potential smoothly follows the corresponding squared one (Figure []]). 
While the continuous potential is defined by five parameters (plus the cutoff distance and 
two additional constants), the squared potential is defined by only three parameters: Ur/Ua, 
the repulsive width Wr/cl = (RR — a)/a of the shoulder and the attractive width wa/o, of the 



U (r) 



Ur 

1 + exp [A (r - R R ) /a] 



— Ua exp 



[r - RaY 
28\ 



3 




FIG. 1: Isotropic pairwise interaction potential in Eq. ([I]), for the parameters in the text (continuous 



50, 



51 



S2| for Ur/Ua = 2, 



line), con spared with the squared potential (dashed line) of Kefs. Q, 
repulsive width wr/cl = 0.6 and attractive width wa/cl = 0.75. The two potentials have the 
same Ur/Ua, the same repulsive range wr = (Rr — a) and the same attractive width wa = 



25a\/2 hi(2). The continuous potential crosses the lines of the shoulder and well of the squared 
potential approximately in the middle. 



well. A straightforward correspondence between the two types of potential can be established 



by choosing the same Ur/Ua, the same wr = (Rr — a) and wa = 25^ a/2 ln(2), respectively. 

The squared potential is more appropriate for a systematic study of the phase diagram 
changes for different combinations of parameters. This analysis has been performed in 



Refs. ]5 



51 



521 ] . The continuous potential, on the other hand, resembles effective pair 



potentials obtained from the second-order pseudo-potential perturbation theory for liquid 
alkaline-earth metals, such as magnesium, calcium, strontium and barium near the melting 



point [55|, [56] , or effective potentia 
by inversion of structural data [57 ]. 



s for colloidal mono-layers at the water-air interface given 
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Moreover, the angle-average of several asymmetric water potentials give rise to isotropic 
soft-core effective interactions, such as the one studied in Refs. 58|, |59j . Hence, the contin- 
uous potential appears to be a more realistic model for different problems. 



III. MOLECULAR DYNAMICS SIMULATIONS 

We perform Molecular Dynamics (MD) simulations in the NVT ensemble for N = 1372 
particles of unit mass m in a cubic box of volume V with periodic boundary conditions 
at temperature T. We use standard MD technique for continuous potentials, with velocity 
Verlet integrator 60J. We keep T constant by rescaling the velocities at each time step 



by a factor of (T/T) 1//2 where T is the current kinetic temperature and T is the desired 



thermodynamic temperature 



60j. We calculate the pressure P by means of its expression 



in terms of the second virial coefficient 60] 



We consider 28 densities p = N/V in the range 0.05 < pa 3 < 0.33 and 5 temperatures in 
the range 0.4 < k B T/U A < 1-9, where k B is the Boltzmann constant. For each state point 
we run 8 independent simulations from different configurations equilibrated at T + 8T, with 
k B 5T/U A = 0.1 for the state points with 0.4 < k B T/U A < 0.8, and with k B 5T/U A = 0.2 for 
the state points with 0.9 < k B T /U A < 1.9. 

Each run has 10 6 MD steps, with a time step 5t = 3.2 x 10~ 3 times units, where time is 
measured in units of (a 2 m/V a) 1 / 2 (of the order of ~ 2.1 x 10 _12 s for for argon-like atoms 
and ~ 1.7 x 10~ 12 s for water-like molecules). We record pressure and potential energy 
every At = 100 MD steps. For the calculation of the averages we consider only the data 
corresponding to states that are stable within the simulation time, as follows from the 
analysis of the time series of pressure, average temperature, and potential energy. Usually 
an equilibration time of 10 5 MD steps is enough to reach a stable state. We average the 
data binned in 10 3 consecutive MD steps, a time that we test to be enough to decorrelate 
the data for the considered state points. 

For densities 0.275 < pa 3 < 0.330 at k B T/U A = 0.5, and for pa 3 > 0.21 at k B T/U A = 0.4, 
we observe spontaneous crystallization during the simulation time. To get data for the 
metastable liquid phase in these cases, instead of annealing from an equilibrium configuration 
at T+5T, we perform instantaneous quenches from configurations equilibrated at k B T/U A = 
0.8. Next, we analyze the evolution in time of the structure factor, as described in details in 



5 



Ref. 



50j, and we use for the averages only the data corresponding to the metastable liquid 



state. We accumulate statistics for these state points up to reach a number of independent 
metastable configurations comparable to the number used for the other state points. We try 
to analyze in the same way also data for quenches to k B T '/V \ = 0.3 and lower T, but the 
crystal nucleation time observed in these cases is too small to allow us reliable estimations. 



IV. THE PHASE DIAGRAM 

We analyze the isotherms in the P-p phase diagram. At low density (Figure. [2]) we 
observe non-monotonic isotherms, showing a maximum at low p and a minimum at higher p 
(the so-called loop), typical of a first-order phase transition between two coexisting phases. 
By joining these minima and maxima we get an estimate of the two branches of the spinodal 
line, i. e. the limit of stability of one of the two coexisting phases with respect to the other. 
The point where the two branches of the spinodal line meet (the vertex of the spinodal) 
is, by definition, the gas-liquid critical point Ci, with critical pressure Pc x , density pc x 
and temperature T Cl . Since, by definition, at pc x the p-derivative of the critical isotherm 
P(p, Tcj) is zero, we interpolate our data to find the isotherm that appears to be flat 
around the approximate vertex of the spinodal line. Our final estimate of the parameters 
for the gas-liquid critical point is p Cl a 3 = 0.08 ± 0.02, P Cl a 3 /U A = 0.021 ± 0.008, and 
k B T Cl /U A = 0.96 ±0.06. 

At higher density we observe that isotherms for k B T '/Ua < 0.8 show a clear minimum 
in the derivative of P with respect to p (Figure [3]) approaching zero at k B T '/U \ = 0.5 and 
pa? 0.25. If we assume that the isotherms below this temperature behave in a regular way 
and with continuity with respect to the isotherms at higher T, then Figure [3] shows that the 
isotherm at k B T '/U \ = 0.5 is just above the critical isotherm, where (dP/dp) T = 0. This 
is confirmed by the analysis of the state points equilibrated at k B T '/U A = 0.4 (small open 
circles in Figure H]) showing a clear maximum in pressure at pa 3 w 0.21, as expected for an 
isotherm entering into the instability region delimited by a spinodal line. Direct observation 
of the loop along the isotherm at k B TU /Ua = 0.4 is hindered by the fast nucleation of the 
crystal, as described in the previous section. 

To estimate the parameters of this high-density critical point C2, we extrapolate from 
the available isotherms at higher T, the high-p state points for k B T/U A = 0.4. The re- 
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suiting curve [the lowest continuous line in Figure HJ^a) and Figure H](b)] is in very good 
agreement with the state points calculated at 0.20 < pa 3 < 0.22, confirming the validity of 
the assumption of continuity and regularity with respect to the isotherms at higher T. 
Therefore, our simulation show a critical point Ci between two liquids at high p and below 



the gas-liquid critica 
studied in Refs. 14 



temperature, as seen for the squared soft-core shouldered potential 



50, |5l|, [52| (Figure [I]). The high density liquid is metastable with 



respect to the spontaneous crystallization and the crystal nucleation process is very fast 
around C<i- The interpolation of the data allows us to estimate the critical parameters of 
the liquid-liquid critical point C 2 as p C2 a 3 = 0.248 ± 0.007, P C2 a 3 /U A = 0.286 ± 0.006, and 
k B T c jU A = 0.49 ±0.01. 



V. COMPARISON WITH THE PHASE DIAGRAM OF THE SQUARED POTEN- 
TIAL 

We compare the phase diagram found here with the phase diagram for the corresponding 
squared potential (Figured]) with parameters Ur/Ua = 2, wr/ci = 0.6 and wa/cl = 0.75. 
The phase diagram of the squared potential with the straightforward correspondence to the 
one studied here is not available in literature, however, in Ref. we studied a squared 
potential with set of parameters Ur/Ua = 2, wr/cl = 0.6 and w^/a = 0.7 (see Figure 3 of 
Ref. [5_1|). Comparison with this potential is shown in Table [H 

TABLE I: Temperatures Tc x an d Tc 2 , pressures Pc x and Pc 2 , and densities pc x and pc 2 , for the 
critical points C\ and C2, respectively, computed for the continuous potential presented in this 
work and the squared potential with parameters Ur/Ua = 2, wr/cl = 0.6 and wa/cl = 0.7 in 
Ref. [51]. 



Potential 


k B T Cl /U A 


a 3 P Cl /U A 


a 3 Pd 


k B T C2 /U A 


a 3 P C2 /U A 


a 3 pc 2 


Continuous 


0.96 ± 0.06 


0.021 ± 0.008 


0.08 ± 0.02 


0.49 ± 0.01 


0.286 ± 0.006 


0.248 ± 0.007 


Squared 


1.24 ±0.01 


0.03 ±0.01 


0.09 ± 0.02 


0.69 ± 0.02 


0.11 ±0.02 


0.28 ±0.02 



The comparison shows that the critical 
are almost the same in the two cases with 
percentage of variation is 11% for pc ± , 23% 



parameters of the gas-liquid critical point C\ 
lower values for the continuous potential. The 
T Cl , and 30% for P Cl . 
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For the liquid-liquid critical point C2 the parameter pc 2 of the continuous potential 
decreases of the same percentage 11% of the critical density of Ci, and Tc 2 decreases of 
29%, a percentage comparable to the decrease of T Cl . The situation is different, instead, for 
the variation of Pq 2 , that increases of 61% in the continuous case with respect to the case 
of the squared potential. 

Therefore, the largest difference between the squared and the continuous case is on the 
pressure of the critical points, with a considerable increase of Pc 2 ■ These effects are consistent 



with the analysis performed in Refs. 5l|, [52j for a decrease of the attractiveness of a squared 



potential, as a consequence of the decrease of the range of the attractive well wa, or the 
increase of the repulsive energy Ur/Ua- Panels (a), (b), (c), (g), (h), (i) of Figure 9 of 



Ref. 



51] , indeed, show that a decrease of wa, or an increase of Ur/Ua, determines the 



increase of Pq 2 and the decrease of Tc 2 , while the other critical parameters for C\ and C2 
are almost unaffected. 

Hence, the continuous potential studied here behaves like a squared potential less attrac- 
tive than the one that would be in straightforward correspondence with it. This is consistent, 
with the observation that the volume integral of the attractive part of the present poten- 
tial is smaller than the volume integral for the corresponding squared potential, while the 
volume integral for the repulsive part is approximately the same in both potentials. 



52j, this continuous potential with two critical points satisfy the 



As predicted in Ref. 
empirical relation 

where Vsc = ^~Rr- With the present choice of parameters we find that the integral in 
Eq.Q has a numerical value of ~ —0.88, supporting the idea that the Eq.(j3J) is a good 
empirical relation to predict if an isotropic potential could display a phase diagram with two 
critical points. 



VI. DENSITY ANOMALY 

A closer look at the high-density part of the P—p phase diagram [Figure H^b)] shows that 
at low T, close to C2, the isotherms cross at various densities. This is a clear signal of a 
density anomaly. 

Indeed, the crossing of isochores implies (dP/dT) v = 0, i.e., for one of the Maxwell rela- 
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tions, (dS/dV) T = 0, where S is the entropy. Hence, is (dS/dP) T = (dS/dV) T (dV/dP) T = 
0, where we use the fact that (dV/dP) T < is finite for liquids. Hence, for one of the Maxwell 
relations, is (dV/dT) p = 0, or equivalently (dp/dT) p = 0. Since at high enough T the iso- 
baric density increases on decreasing T, this relation implies that there is a temperature of 
maximum density Tmd at constant P (Figure ED- 

For T > Tmd the isobaric density decreases for increasing T, as for normal liquids, while 
for T < Tmd it decreases for decreasing T, giving rise to an anomalous behavior in density. 
This can be seen in a clear way by studying directly the (dP/dT) p (Figure [6]). Comparison 
with Figure H](b) shows that the isochores with a minimum, i.e. (dP/dT) p = 0, correspond 
to the densities where the isotherms cross in the P-p phase diagram. 



VII. DISCUSSION AND CONCLUSIONS 



Our results show that the continuous interaction potential presented here has a maximum 
in density along the Tmd line. This anomaly in density is typical of network-forming liquids, 

J6l|, and has been observed in other soft-core potentials 



such as water or SiO 



39 



40 



41 



43 



44 



45 



46. 



43, 



481 ] . In the present model, as well as in those for water and 



other anomalous liquids 29] the Tmd hne converges toward a liquid-liquid critical point, i.e. 
C*2 in this work. This observation has been often considered as an evidence for the fact that 
the liquid-liquid critical point implies the presence of a Tmd line. 

On the other hand, we know that for the squared soft-core potential the liquid-liquid 
critical point occurs in absence of the density anomaly, as shown by the extensive analysis 
performed in Ref. 50]. As it was observed (see Figure 21 in in Ref. jsO]), the squared 
potential does not satisfy the condition for the density anomaly, Eq.(21) in Ref. 50], but 
approaches that condition in an asymptotic way. Hence, is not clear if the absence of the 
Tmd is a consequence of the specific discontinuous shape of the squared potential. 

Since the continuous potential presented here can be tuned to approximate in a close way 
the discontinuity at Rr for the squared potential (Figure [7]), it could be interesting to study 
the effect on Tmd of the variation of the parameters A (determining the slope at the Rr) 
of the continuous potential. Such a study could shed light on the relevant question about 
the relation between the liquid-liquid critical point and the presence of the density anomaly. 
It is worth to observe here that the implication of the presence of a liquid-liquid critical 
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point for network-forming molecular liquids with anomalous density has been supported by 
extensive studies on water-like models (see for example Refs. 16, [l?], Q)- 

In conclusion, we have presented here a continuous soft-core potential with a shape which 
resembles that of effective interaction potentials for liquid metals [3] or colloid solutions 



57| . The model displays a phase diagram with a high-T gas-liquid critical point C\ at low 



49 



50 



51 



52| shows that 



density, and a low-T liquid-liquid critical point C2 at high density, occurring at the limit of 
stability, within the simulation time, of the liquid phase with res pec t to the crystal. 

Comparison with the case of the soft-core squared potentials 
the continuous potential is less attractive than its straightforward corresponding squared 
potential. The main qualitative difference with the case of squared potentials is that the 
continuous potential shows the density anomaly, typical of network-forming liquids such as 



39, 4(j, 



41 



43 



45 



46 



47 



481 ] . Further 



water, and observed in other soft-core potentials 
investigation is necessary to understand how the appearance of this anomaly is related to 
the parameters of the potential. 
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FIG. 2: Low density phase diagram with the estimate of the gas-liquid critical point C\ (large full 
circle) at p Cl a 3 = 0.08 ± 0.02, P Cl a 3 /U A = 0.021 ± 0.008, and k B T Cl /U A = 0.96 ± 0.06. Symbols 
are state points calculated with MD in the pressure-density (P-p) plane for, from top to bottom, 
k B T/U A = 1.9, 1.7, 1.5, 1.3, 1.1, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4. Full lines are cubic interpolations of 
the isotherms. Dashed lines are a guide for the eyes connecting symbols for the isotherms at which 
the cubic interpolation is not reliable. Dotted lines are liner interpolations between the isotherm 
at k B T/U A = 1.1 and the isotherm at k B T/U A = 0.9 for (upper thin dotted line) k B T/U A = 1.02 
and (lower thick dotted line) k B T/U A = 0.96, approximating the critical isotherm. The thick gray 
dashed line at T < Tc 1 is a guide for the eyes connecting the points of minimum P along the 
isotherms for p > pc x (i. e. the limit of stability of the liquid with respect to the gas), and the 
points of maximum P for p < pc x (i- e. the limit of stability of the gas with respect to the liquid), 
and representing an estimate of the spinodal line. Errors on state points are smaller then symbols 
size. 
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FIG. 3: Numerical derivative of P with respect to p along the isotherms at, from top to bottom, 
kftT jU a = 1-9, 1.7, 0.8, 0.7, 0.6, 0.5. Symbols are as in Figure[2j At high T the derivative increases 
with p as for normal liquids, while at low T the derivative is non monotonic with a clear minimum 
that approaches zero for ksT/UA = 0.5 at density pa 3 « 0.25. Errors on state points are smaller 
then symbols size. 
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FIG. 4: (a) Complete phase diagram with the gas-liquid critical point C (large full circle) as 
shown in Figure [2j and the liquid-liquid critical point C2 (large open circle) estimated at pc 2 a 3 = 
0.248 ± 0.007, P C2 a 3 /U A = 0.286 ± 0.006, and k B T C2 /U A = 0.49 ± 0.01. Symbols and lines are 
as in Figure [2j The lowest solid line below C2 is an extrapolation of the data at ksT/UA = 0.4 
in the region of inevitable crystallization, on the basis of the interpolation of the isotherms at 
higher temperatures. This extrapolation shows very good agreement with the state points that 
can be directly calculated at k^T jXJ a = 0.4 (small open circles). In the high-p region, from top 
to bottom, the upper thin dotted line is the interpolation at k^T j\J a = 0.55 between consecutive 
isotherms, the thick dotted line is the interpolated critical isotherm, and the lowest thin dotted 
line is the interpolations for fc^T '/V 'a = 0.45. The thick gray dashed line at T < Tc 2 is a guide for 
the eyes connecting the points of minima and maxima of P along the isotherms representing an 
estimate of the liquid-liquid spinodal line, (b) The high-p portion of the phase diagram, around 
the liquid-liquid critical point C2. Symbols and lines are as in the previous panel. For sake of 
clarity we show only the lines corresponding at the interpolations of the data along the isotherms 
for 0.4 < ksT/UA < 0.8. The agreement between the data at k B T '/U A = 0.4 (open circles) and the 
extrapolation of the corresponding isotherm (lowest solid line) on the basis of the data at higher 
T is remarkable. In both panels, when not shown, errors are smaller then symbols size. 
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FIG. 5: The P—T phase diagram shows isochores with a clear minimum. The gas-liquid critical 
point C\ and the liquid-liquid critical point C2 are shown on the basis of the estimates in Figured! 
They are consistent with the crossing occurring at the highest T along each isochore. Symbols 
correspond to isochores for, from bottom to top at ksT/UA = 1, pa 3 = 0.005, 0.075, 0.100, 0.125, 
0.150, 0.175, 0.200, 0.205, 0.210, 0.215, 0.220, 0.225, 0.230, 0.235, 0.240, 0.245, 0.250, 0.255, 0.260, 
0.265, 0.280, 0.290, 0.300, 0.310. Dashed lines are a guide for the eyes. Solid lines are a quartic 
polynomial fit between the isochores showing a minimum. The open circles represent the point 
where (dP/dT) p = for the non-monotonic isochores. The thick gray line is a guide for the eyes 
representing an estimate of the Tmd temperature of maximum density at constant P. When not 
shown, errors are smaller then symbols size. 
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FIG. 6: Derivative of pressure with respect to temperature along the isochores at (from bottom to 
top) p = 0.200, 0.205, 0.210, 0.215, 0.220. Lines are plots of the cubic functions resulting from the 
derivative of the best-fit polynomials in Figure EJ with a relative errors on the parameters smaller 
than 10~ 3 . The temperature T* where the derivative is zero is the minimum along the isochore, 
with T* = Tmd at the pressure corresponding to the state point at the given density and T*. 
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FIG. 7: The continuous potential can approximate the discontinuity at the repulsive distance 
Rr/cl = 1.6 of the squared potential (thin dashed line) for increasing values of the parameter A, 
determining the slope at Rr. Continuous line is for A = 15, studied in this work. Dash-dotted 
line is for A = 30. Thick dashed line is for A = 100. 
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